#ifndef RTNORMNC_HPP
#define RTNORMNC_HPP


//Implementation of Robert's truncated Gaussian (1995)
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include "univariate/Main/ran_tg_zig1000.h"
#include "../Random.hpp"
#define pINF 10000

class rtnormNC
{
	public:
	rtnormNC()
	{
		int seed = RandomG::getSeed();    
		//T = gsl_rng_ranlxd1;  // generator with good math properties
		//T = gsl_rng_taus; //default generator
		T = gsl_rng_mrg;
		r = gsl_rng_alloc(T);    
		gsl_rng_set (r, seed); 	  
		boost::random::uniform_01<> U;
		unif=new RandomG::Random<boost::random::uniform_01<> >(U);
	}
	~rtnormNC()
	{
		//delete unif;
		delete T;
		delete r;
	}
	double operator()(double m,double mm,double mp,double s)
	{
		double res=0;
		if(-mm==pINF)
		{
			
			mp=(mp-m)/s;			
			mm=-mp;
			double x = ran_tg_zig1000(r,mm);	
			res=(-x*s+m);	
		}else if(mp==pINF){
			int i=0;
			mm=(mm-m)/s;
			double x = ran_tg_zig1000(r,mm);	
			res=(x*s+m);	
		}else{
			double i=0;
			while(i==0)
			{
				double z=(mp-mm)*(*unif)()+mm;
				double r=rho(z,mm,mp);	
				double u=(*unif)();
				if(u<=r)
				{
					i=1;
					res=z*s+m;				
				}
			}
		}
		return res;
	}
	double rho(double z, double mm, double mp)
	{
		double res=0;
		if(0>mp)
		{
			res=exp((pow(mp,2)-pow(z,2))/2);
		}else if(mm>0){
			res=exp((pow(mm,2)-pow(z,2))/2);
		}else{
			res=exp(-(double)pow(z,2)/2);
		}
		return res;
	}

	private:
	RandomG::Random<boost::random::uniform_01<> > *unif;
	const gsl_rng_type *T;
	gsl_rng *r;
};

#endif

